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I. STATEMENT OF THE PROBLEM 


1.1 Problem Definition and Method of Attack 

The objective of this research work is to numerically simulate the vacuum plume 
flow field in the backflow region of a low thrust nozzle exit. In space applications, the 
low thrust nozzles are used as a propulsion device to control the vehicle attitude, or to 
maneuver the vehicle flight trajectory. When the spacecraft is deployed in the orbit or 
cruising in a planetary mission, the vacuum plume is created behind the nozzle exit (so 
called backflow region), by the exhausting gas of the propulsion system or by venting 
internal gas to the extremely low density ambient. The low density vacuum plume flow 
regions cover the continuum, transitional and free molecular flow regimes (see Figure 
1), which were characterized by the Knudsen number K n , 

K — ' Km 

hn ~T 

where A m is the mean free path of the gas molecules and L is the characteristic length 
of the flow field. Figure 1 shows the backflow regions determined by the A n number. 
The transitional regime is defined by 0.01 < K n < 10. The conventional Navier-Stokes 
equations are valid only in the flow region close to the nozzle exit since the validity of 
the Navier-Stokes equations fails asymptotically as the Knudsen number increases. The 
vacuum plume characteristics prediction is primarily a problem of transitional aerody- 
namics. 

In order to handle the rarefaction effects properly, the following hypothesis is made 
in this study, 

“the Burnett equations are the governing equations in the tran- 
sitional rarefied aerodynamic regime, where the Knudsen num- 
ber is approximately between 0.01 and 10.' 1 
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This hypothesis is a compromise between the Boltzmann equation and the Navier Stokes 
equations. Based on this hypothesis, the method to attack the vacuum plume problem in 
the backflow region is to develop a numerical technique to analyze the characteristics of 
the vacuum plume by solving the Burnett equations. For clarity, it is further assumed 
that the Navier-Stokes equations are no longer valid for this research problem, even 
though the Navier-Stokes equations are still approximately correct up to K n = 10 
(Cheng, 1989). One of research objectives is to verify the consistency of the Burnett 
solutions with Navier-Stokes solutions and the Direct Simulation Monte Carlo solutions 
as well as the source flow correlation method in the transitional regime. 
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Figure 1. Schematic flow regimes in the backflow region and Validity of the Algorithms. 
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1.2 Survey of Low Density Nozzle/Space Plumes 

Over the past years, many research efforts have been put into the space plume studies 
because the gasdynamics of this rarefied plume affects the vehicle trajectory control 
and propulsion efficiency. The plume in backflow region creates a negative force to the 
space vehicle, affects the vehicle torque control (Baerwald, 1978). Also the exhausted 
plume may contaminate the surface of the spacecraft. For the long duration flight, such 
as the Moon/Mars mission, an accurate definition of the vacuum plume is required for 
the improvement of the spacecraft performance. 

It is difficult to simulate the space plume conditions by the ground test facility since 
it is very hard to pump the testing vacuum chamber pressure lower than 1 micron. It 
has been known that the performance of the low-thrust nozzle, unlike high thrust nozzle, 
is greatly affected by the back pressure in the testing vacuum chamber, even though the 
nozzle is operating under “choked” conditions. Available wind tunnel data show that, 
for a micropound thruster, if the back pressure falls into the region of 10 to 1000 microns, 
the thruster performance decreases with decreasing back pressure, however, the further 
decreasing of back pressure in the range of 0.1 to 10 microns will increase the thruster 
performance, (Pugmire 1968, John 1966, Sutherland 1966). There has been no proven 
theory in the exact cause of this “low pressure” anomaly. The numerical prediction is 
the only way to provide the answer. 

The low density plume possesses the following characteristics (Kogan 1986, Legge 
1988, Bird 1976): (1). The flow region covers the continuum, transitional, and free 
molecular flow regimes without solid boundaries; (2). The conventional Navier-Stokes 
equations are only valid close to the nozzle exit; (3). The translational nonequilibrium 
effects occur in the transitional regime; (4). There is no flow separation; (5). Transport 
properties not only depend on the first order spatial gradients of macroscopic flow 
variables, but also high order spatial gradients. These information indicate that the 
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vacuum (space) plume prediction is a problem of the transitional aerodynamics. 

There have been many efforts towards solving backflow vacuum plume problems. 
The source flow’ technique developed by Boynton (1967) and Simons (1972) is a useful, 
inexpensive engineering design tool for the calculation of spacecraft impingement effects 
in the continuum regime, although some characteristics of the expansion are not properly 
treated. In order to produce a realistic source flow plume, the maximum Prandtl- Meyer 
expansion angle 0 max is computed at the nozzle lip. Using this angle, the source flow 
area is computed as assuming the source to be located on the axis at a point defined by 
the interaction of the expansion angle with the axis. The model primarily assume that 
in the backflow plume region, the density at some radial distance from the source may 
be expressed as 

P_ _ ± 7 ( * 9 \ 

p 0 r 2 C ° S \2 0 m ax) 

where A is a plume constant, and, velocity may be assumed constant. The exponent 7 
is obtained from the conservational mass and momentum. At the large distance from 
the nozzle exit, the flow may fall into transitional and free molecular regime, the gas is 
no longer in thermal equilibrium and the continuum equations become invalid. In order 
to understand detailed fluid dynamic phenomena in the backflow region, this empirical 
method can not produce successful results. 

Method of Characteristics (MOC), plus the boundary layer prediction technique, 
can predict accurate plume expansion if Mach number is less than 10. For vacuum 
plume, the density is very low, the conventional MOC technique can still predict the 
flows, but the physical meaning and accuracy out of the computational number is in 
doubt. 

Navier- Stokes equations have been widely accepted as a powerful and accurate 
governing equations for most continuum fluid flows. Most of the Navier-Stokes plume 
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predictions are performed at the low altitude where the density is fairly large compared 
to vacuum plume. Vacuum plume prediction is a problem of highly nonequilibrium 
transitional aerodynamic problem. In the discussion later, the Navier-Stokes equations 
fails in the transitional regime. 

It is necessary to have new computational concepts (Direct Simulation of Monte 
Carlo technique) or governing equations (such as Burnett equations with nonequilibrium 
effects) in order to solve the complex transitional aerodynamics (vacuum plume in the 
backflow region), or a new numerical technique to solve the difficult physical equations 
(Boltzmann equation). 
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II. TECHNICAL DISCUSSION 


2.1 Introduction 

From the physics point of view, the Boltzmann equation has long been established 
as the standard mathematical formulation of a nonequilibrium thermally perfect gas as a 
set of molecules. Whereas the continuum Navier-Stokes equations have the flow velocity 
and macroscopic thermodynamic properties as dependent variables. The only dependent 
variable in the Boltzmann equation is the distribution function for the molecular states. 
The velocity distribution function / is a scalar function which depends in general on 
time, the three components of molecular velocity, and the three components of molecular 
position, 

f(u,,x,,t)du t = probability of molecule in velocity range Uj + du,, 

at position r;, at time t 

By taking moments (Lumpkin, 1990) of the velocity distribution function over all ve- 
locity space, all macroscopic thermodynamic properties of the gas can be obtained. 
However, the number of dependent variables is reduced at the expense of the additional 
number of independent variables from those of physical space to those of mathematical 
phase space. This leads to difficulties in the way of direct numerical or analytical solu- 
tions of Boltzmann equation for nontrivial gas flow problems though analytical solutions 
exist for a few simple cases. This applies even for monatomic gas flows and there is no 
prospect of obtaining direct solutions of the Boltzmann equation for complex aerother- 
modynamic problems involving real gas effects. 

One alternative is to model the gas flow at the molecular level and a number of 
Monte Carlo simulation methods have been developed to do this. The Direct Simu- 
lation Monte Carlo (DSMC) method has gained a broad acceptance. Yet, while this 
method is capable of good predictions of hypersonic flow in the continuum transitional 
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regime, it is quite expensive in terms of computer time because the computational time 
is proportional to the molecules and cells number in the simulated flow. DSMC assumes 
that molecules are point centers of the mass and point centers of repulsion containing 
no internal structure. DSMC is primarily employed in highly rarefied gas flows, it is not 
practical to use DSMC for high Reynolds number flows in which Navier-Stokes equations 
can be solved instead. As a numerical technique, DSMC may capable of predicting com- 
plex problems involving effects such as dissociation, chemical reactions, and radiation 
as long as the physics of the various phenomena can be described at the molecular level. 
However some essential physical questions concerning the continuum and transitional 
aerodynamic phenomena such as entropy are not answered satisfactorily. 

To circumvent the vexed Boltzmann equation, an alternative can be achieved by 
solving conservational equations. To relate the nonequilibrium kinetic theory (Boltz- 
mann equation) to the continuum theory of gas dynamics, Hilbert-Chapman-Enskog 
expansion method plays an important rule. Based on this method, the velocity distri- 
bution function can be expressed in the series expansion of Knudsen number, 

/ =/(«> + /<» + /< 2 > + + ■ • ■ 

= fi 0) aiK n + 02-^n "t" a 3 ~t~ ’ ' ') i 

which is the perturbation expansion of the velocity distribution function about the 
Maxwellian distribution which a gas exhibits at equilibrium, where is the Maxwellian 
distribution function, a t , i = 1,2, 3, ... are functions of density, molecular velocity, and 
temperature (Bird, 1976). Substituting this expansion into Boltzmann equation, taking 
moments of the Boltzmann equation yields the continuum equations of fluid mechanics, 
a set of conservation equations describing global conservation of density, momentum, 
and energy. To close this system of equations requires constitutive equations which 
express viscous stress and heat flux in terms of the distribution function rather than 
macroscopic gradients. Closure of these equations can be achieved by using any of 
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the nth order approximations to the distribution function described in the expansion. 
The constitutive conditions relating viscous stress tensor and heat flux to macroscopic 
gradients can be obtained. By equating like terms for Knudsen number in Boltzmann 
equation, it was proved that the Euler equations are the zeroth-order approximation, 
the Navier-Stokes equations are the first-order approximation and the Burnett equations 
are the second-order approximation. 

2.2 Inadequacy of the Navier-Stokes Equations 

The conservation equations of fluid mechanics are valid in all flow regimes. But 
solving the general conservation equations (such as Navier-Stokes equations) require 
constitutive relations to close the system of equations. When the local Knudsen number 
which is based on the scale length of the gradients of the macroscopic flow properties 
exceeds 0.1, the conventional constitutive equations that relate the shear stress and heat 
fluxes to these gradients breakdown. One of the major breakdown in the transitional 
regime is that the thermal equilibrium fails. 

Reasons for the failure of the Navier-Stokes equations may be associated with any 
or all of the following assumptions embodied in these equations (Fisko, 1989): 

(1) . Linear stress-strain tensor dependent only on velocity gradients. 

(2) . Linear heat -flux vector dependent only on temperature gradients. 

(3) . Zero bulk viscosity. 

(4) . Small Knudsen number flow (Continuum gas model). 

(5) . No direct contribution of nonequilibrium internal molecular energy to viscous stress 

or heat flux. 

2.3 Difficulties of Solving Burnett Equations 

It has been a mislead concept that the Burnett equations are no more accurate than 
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the Navier-Stokes equations until Fisko and Chapman (1988) obtained the Burnett so- 
lutions for thick shocks. The comparison of solutions of Navier-Stokes, Burnett and 
the DSMC clearly shows that the Burnett solutions are superior to the Navier-Stokes 
solution (Figure 2). However, solving the Burnett equations is by no means a straight- 
forward procedure because there are not enough well-defined boundary conditions for 
the fourth-order nonlinear partial differential equations. In addition, the questions of 
convergence, accuracy, stability, and viscous dissipation, etc., still remain. 

Several fundamental questions are addressed in this report: 

1. For tested monatomic gases and tested Mach numbers, the Burnett equations give a 
significant improvement over the Navier-Stokes equations in all macroscopic shock 
parameters as compared with DSMC and available experimental results (Lumpkin, 
1990). But the flows in the transitional regime are always nonequilibrium. What is 
the theoretical ground to accept that the Burnett equations are the true governing 
equations in the transitional regime? The same question was raised towards the 
Boltzmann equation. 

2. The mathematical characteristics of the Burnett equations are not clear, primarily 
because they are third-order highly nonlinear partial differential equations. 

3. Most of the numerical schemes were developed for the Euler and Navier-Stokes 
equations. Should we routinely apply the numerical scheme to discretize the Bur- 
nett equations? What is the best stable numerical scheme to solve Burnett equa- 
tions efficiently? 

4. The boundary conditions for solving Burnett equations are undetermined. 

As for this report, we are still working on the numerical scheme formulations and 
one-dimensional shock solutions in the wind tunnel conditions without answering these 
basic questions. 

As a numerical experiment on Burnett equations, Direct Lower-Upper Factorization 
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(DLUF) scheme and the Arbitrary-Lagrangian-Eulerian (ALE) scheme have been devel- 
oped and adopted to obtain the solutions for Burnett equations. The one-dimensional 
normal shock wave structures have been obtained by solving Navier-Stokes and Bur- 
nett equations, unsteady shock-tube solutions are also obtained from the ALE scheme. 
Difference between Navier-Stokes and Burnett solutions are shown for the density- 
temperature profiles for normal shock structure. The developed DLUF scheme shows 
a good potentials for solving Burnett equations will success. The ALE scheme has 
the advantages of reducing the order of partial differential equations, and allows more 
flexibility for viscous stress tensor and heat flux formulations. 



Figure 2. Comparison of reciprocal shock wave thickness by Navier-Stokes, 
Burnett, DSMC solutions and experiment. 
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2.4 Translational thermodynamic Nonequilibrium Effects 

One of the major characteristic of vacuum plume flows is that the flow is always nonequi- 
librium, which primarily include the molecular translational, vibrational, and rotational 
nonequilibrium. For monatomic gases, the rotational and vibrational nonequilibrium are 
absent. Burnett equations do not account for molecular internal energy (i.e. rotation or 
vibration), and thus it is applied for monatomic gases. Since most of the gas is diatomic, 
attention must be given to the theory of the rotational and vibrational nonequilibrium 
effects. Unfortunately, the theory for rotational nonequilibrium is not as mature as the 
theory for translational nonequilibrium. This is due to the facts of complexity of the 
problem and lack of study interest. 

In view of the theories up to date, it would seem that there are two basic ways 

to model rotational nonequilibrium effects using a continuum formulation. The first 

* 

method, the so called “bulk viscosity” method, uses modifications to the stress ten- 
sor to account for the dissipation which results from rotational nonequilibrium. It 
assumes rapid equilibration and small departures from equilibrium, and is therefore in- 
appropriate for hypersonic flow conditions. The second method, so called “relaxation 
equation” method, uses two temperatures to describe the rotational and translational 
energy modes. Since an additional temperature is used to describe the macroscopic 
state of the fluid, an additional governing equation is needed in addition to the conser- 
vation equations for mass, momentum, and energy. The additional governing equation 
can be derived from either the fully classical or the quasi-classical viewpoints. In both 
cases, the assumption of long relaxation times is required. The quasi-classical derivation 
also requires the assumption of small departure from equilibrium, an assumption which 
is not strictly valid in hypersonic flow. Comparing these two methods, for low den- 
sity hypersonic flows, the relaxation equation method is preferred over a bulk viscosity 
method. 
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The general theory of translational-rotational nonequilibrium has only been car- 
ried out for first order departures from equilibrium distributions and in the limit of long 
relaxation times. The resulting set of equations (use both classical and quasi-classical 
models) is simply the Navier-Stokes equations supplemented by an additional govern- 
ing equation, such as Jean’s equation (Jean, 1904). The implementation of rotational 
nonequilibrium into Burnett equations, even for a monatomic gas, with the simplified 
modifications, is uncertain because of the accuracy of the additional governing equations 
is unknown for low density hypersonic flows. 

the rotational and vibrational thermodynamic nonequilibrium excited in the high 
temperature, hypersonic flows are not important. 

In this study, the interest is in the simulation of the backflow region, where low 
density effects is the major governing mechanism, the rotational and vibrational ther- 
modynamic nonequilibrium excited in the high temperature hypersonic flows are not 
important. Therefore, our intention is to solve the monatomic gas problem such that 
translational nonequilibrium is the most important phenomena required to be accurately 
implemented. There are several attempts to analyze the translational nonequilibrium 
effects (Anderson, 1989). To handle the translational nonequilibrium, the translational 
and kinetic temperature have been introduced. Similarly, non-isotropic “pressure” (nor- 
mal and tangential pressure) are used to evaluate the stress tensor. The key questions 
are still remaining: what is the relationship between these properties? What are the 
physical laws that guarantee the essential validity of these concepts? What are the 
appropriate boundary conditions? 

2.5 Slip- Wall Conditions in Nozzle/Orifice Exit 

Due to the complexity of the physical problems, currently, the flows inside the 
low-thrust nozzle and vacuum plume are decoupled, in the sense that, the nozzle exit 
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conditions axe used as the inflow boundary conditions for vacuum plume predictions. 
The specification of the nozzle exit conditions affects the accuracy of the vacuum plume 
predictions. 

Theoretically, boundary conditions and, for unsteady flows, initial conditions must 
be given so as, in principle, to ensure the existence and uniqueness of the solution 
of the problem. However, in practical applications, the existence and uniqueness of 
the solution can not be proved and one must restore to physical intuition, non rigor- 
ous mathematical interpretation, and numerical experiments in order to determine the 
proper boundary conditions. The conventional wall boundary conditions, such as no-slip 
and no-temperature jump conditions, are only valid when K n is less than 0.01. For low 
density nozzle, the nozzle expansion ratios are usually big enough that most likely, the 
flows at nozzle exit have very low densities such that K n is greater than 0.01. In this 
case, the so-called “slip- wall” conditions need to be used at the nozzle wall exit. Based 
on kinetic theory, depending on the thermodynamic state and the laws of interactions 
between the gas molecules and the wall, the slip-wall conditions can be determined 
(Hollanders, 1988). 

The flow variables along body (nozzle wall exit) surface (s) can be computed by 
the first-order Maxwell/Smoluchowski (Kennard 1938, Schaaf 1961) slip boundary con- 
ditions (in Cartesian coordinates), 

Surface velocity 

2 — Ra n 2fi I 7 r / du\ 3 fi / dT\ 

a n p \ 8RT \dy ) s + 4 pT \dx ) s 

and surface temperature T 5 , 

rji _ T m(?L\ 

• w+ a, T+lpP r \SRT\dy), 

where T w is the wall surface temperature, cr n ,a t are the surface reflection and accom- 
modation factors, respectively. For a complete accommodation, cr n = 1, <r t = 1. 
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In order to treat rotational nonequilibrium effects, the above slip-wall boundary 
conditions can be extended as (Kennard 1938), 

Surface velocity u 5 , 


2 — Ra n 2/i I 71 r / du \ 3 p / dT \ 

J\8RT\d^) s + l^T \ fa) s 


surface temperature T 5 , 



3A f 

2pcc vt 



and rotational surface temperature T r5 , 


T r3 — T w + 


2 2A r 

(7^ pCCyr 



where T r is the rotational wall surface temperature, A*,A r are translational and rota- 
tional heat conductivity, respectively, C v t,C vr are translational and rotational specific 
heat per unit volume, and 

c = 



14 - 



III. TECHNICAL ACHIEVEMENTS 


3.1 Reviews of the Burnett Equations 

Literature survey on Burnett equations shows that only limited papers have been 
published on the derivation works. Chapman and Cowling (1970) derived the two- 
dimensional Burnett equations in Cartesian coordinates. Stanford’s group, Chapmann 
and Zhong et al. (1988, 1991), has systematically published their results in solving the 
Burnett equations. 

In order to check the validity of the formulations, Burnett equations have been 
derived for one- and two-dimensional Cartesian coordinates. The final forms of the 
derivations (Liaw, Guo, 1992) is shown in Appendix. 

3.2 Pressure-Based Methodology (ALE Scheme) 

In order to select a best numerical scheme to solve Burnett equations, extensive 
reviews on current existing TVD family of second-order central and upwind difference 
schemes have been performed. It is impossible to adopt a numerical scheme to solve 
Burnett equations without testing the capabilities of the specified scheme. Two different 
kinds of schemes have been examined/developed carefully, pressure-based methodology 
and density-based methodology. 

The Arbitrary-Lagrangian-Eulerian (ALE) numerical technique is adopted in this 
investigation. The ALE scheme is a pressure-based methodology, which is a time- 
dependent finite volume differencing scheme for arbitrarily shaped geometry. In the 
current formulation, the gas could be chosen to be monatomic, diatomic or polyatomic 
arbitrarily. Due to the testing and validation purpose, the high order Burnett stress and 
heat flux terms are switched on only in one direction. For the present time, due to the 
limitation of time duration, the rotational and vibrational nonequilibrium for diatomic 


- 15 - 



ot polyatomic gases are neglected. 

The integrated mass, momentum, and energy equations over a moving control 
volume for gas/fluid may be written as 
Continuity equation: 


d_ 

dt 



pdV + 



u — U) ■ dA = 0 


Momentum equations: 


d 

dt 


/// pudV + f f 

JJJV(t) JJA(t) 


pu{u — U) ■ dA = 



PdA + 



r • cL4 


Energy equation: 


d_ 

dt 


fff pEdV+ff pE(u — U) ■ d~A = — f [ 
JJJv(t) JJA(t) JJcs 

+ fff t : VudV - ff 
JJJv(t) J J A{t) 


Pu ■ dA 


J ■ dA 


where 


2 

J = -kWT - P DJ2 ^V(^) + J<*2) 

m=l ^ 


r = /i( Vu + (Vu) T ) + AV • ul + r (2) 


where p is local fluid density, u is local fluid velocity vector, U is the moving boundary 
velocity, E is the specific internal energy not including chemical energy, P is static 
pressure, r is the viscous stress tensor, r^ 2) is the Burnett stress tensor, J is the heat 
flux vector, P 2 ) is the Burnett heat flux vector, k is the thermal conductivity, A and p, 
are the dilatation and shear viscosity coefficients, respectively. The ideal gas equation 
of state is also assumed for the freestream gas. 
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At each cycle (time step), the calculations is divided into two stages, a Lagrangian 
stage and a rezone- Eulerian stage. In the Lagrangian stage, the velocities move with 
the fluid particle velocity, implicit differencing is used for all the diffusion terms and 
the terms associated with the pressure wave propagation. There is no convection across 
cell boundaries. The coupled implicit equations are solved by the method similar to 
the SIMPLE (Patankar, 1980) algorithm, with individual equations being solved by the 
conjugate gradient method (O’Rourke, 1986). In the rezone Eulerian stage, explicit 
methods are used to calculate convection. The convection calculation can be subcycled 
an arbitrary number of times, and thus the main computational timestep is not re- 
stricted by the Courant stability condition of explicit methods. A second-order upwind 
differencing technique is used for the spatial differencing. Detailed information about 
the ALE scheme can be found in Amsdan (1989) and Deng (1991). 

The implementation of Burnett viscous stress and heat flux terms require special 
treatment of third-order partial differential terms and cross product terms of first-or 
second-order differentials. Due to the integration methodology, the order of partial 
differential terms in the formulations can be reduced by one. It is only necessary to 
define the viscous stress tensor and heat flux vector itself, which include only second- 
order partial differentials and cross product of first-order partial differentials. The basic 
methodology used here is based on the Gauss theorem, for any arbitrary scalar quantity 

Q, 


IIIv dx t dV ~ JJ A Qe ' ' dA ~^ Q ° ei ‘ dAa 

By defining scalar property at the cell center, and vector property at cell vertex, it is 
possible to determine first-order partial derivatives at the cell surface, and second-order 
derivatives at the cell center. 
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3.2.1 Unsteady Shock- Tube Solutions by the ALE Scheme 

To test the validity of the current scheme, a typical one-dimensional shock-tube 
solutions was obtained. The selected gas in this calculations was arbitrarily chosen to 
be Nitrogen. Rotational and vibrational nonequilibrium effects were neglected at the 
present time. The pressure ratio is 10 at the diaphragm initially. The shock tube length 
is lm, the grid point along the tube is 201. The computed shock wave propagation 
at a given time from Euler, Navier-Stokes, and Burnett equations were compared to 
the exact Riemann solutions. The non-dimensional plots of pressure, density, velocity, 
and temperature are shown in Figure 3, 4 and 5. The computed results were based on 
real gas properties, which is currently defined in the way that all the thermodynamic 
properties such as viscosity, thermal conductivity, heat capacity, specific heat ratio, etc. 
are function of local temperature. Results shown at Figure 3 and 4 indicates that the 
current numerical scheme (ALE) predicts the shock wave propagation very accurately. 
The Burnett solutions is very close to the Navier-Stokes solutions in the way that only 
a small difference can be detected from data files. These results validate the current 
formulations. 
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Euler vs. Exact Solution- Shock-Tube Euler vs Exact Solution Shock-Tube 



Figure 3. Comparison of Euler shock-tube solutions with exact soluti 




Navier Stokes vs. Exact Solution Shock-Tube Navier Stokes vs. Exact Solution Shock-Tube 
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X-Positions X-Positions 

* 

Figure 4. Comparison of Navier-Stokes shock-tube solutions with exact solut 






Burnett vs. Exact Solution Shock-Tube Burnett vs. Exact Solution Shock- Tube 



X- Positions . X Positions 

Figure 5. Comparison of Burnett shock-tube solutions with exact solutions. 




3.2.2 Shock Wave Structures by the ALE Scheme 


High Mach number normal shock wave structure prediction is in progress. The 
comparison of Euler, Navier-Stokes, and Burnett solutions will be completed in the near 
future. The numerical experiments are for the Argon (monatomic) gas with different 
high inflow Mach numbers. Different scale calculations will be investigated. 
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IV. FUTURE WORK 


As mentioned in section one, the objective of this investigation is to find the flowfield 
in the vacuum plume backflow region in space. Three interrelated numerical techniques 
are proposed for the upcoming investigation. 

1. The Burnett equations bridge the gaps between continuum and Boltzmann equa- 
tion. The steady state solution of normal shock structure for Argon has been 
obtained by integrating the one- dimensional Burnett equations. The Direct Lower- 
Upper Factorization scheme is used to discretize the governing equations. More 
test conditions, primarily with lower densities and higher Mach numbers, will be 
exercised in this code. Preliminary results show that the Burnett terms suppress 

shock oscillations downstream, and it converges to the steady solution faster. The 

* 

code development for two-dimensional flows will be continued in the future. The 
Direct Lower-Upper Factorization Scheme (DLUF) in two-dimensional rectilinear 
coordinates for Navier-Stokes/Burnett equations with a general grid generator will 
be developed. Under subcontract, Dr. J. D. Mo of the Memphis State Univer- 
sity will be responsible for the basic code development. The principal investigator. 
Dr. G. S. Liaw of Alabama A&M University will be responsible for the insertion 
of the rarefaction effects and the nonequilibrium effects into the basic code. This 
new computer code will calculate the low density space-plume flowfield where the 
Knudsen number is less than 1. The ALE technique is also applied to solve Burnett 
equations. The numerical formulations are in the forms of three-dimensional fluid 
flows for real gases (monatomic, diatomic or polyatomic). Only one-dimensional 
cases has been validated. Continuous code development and validations for two- 
and three-dimensional cases will be performed in the future. The key questions 
here are how to handle the nonequilibrium effects. For monatomic gases, the rota- 
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tional and vibrational nonequilibrium are absent. Most of the real gases composed 
of the vacuum plume are diatomic or polyatomic. The future work towards solving 
Burnett equations including real gas effects will concentrate on the handling of the 
translational, rotational nonequilibrium effects. 

2. For far field where the Knudsen number is larger than 1, the transitional nonequi- 
librium effects dominate the flow characteristics and continuum assumptions are in 
doubt. There are two approaches to tackle this problem: 

(a) . The Direct Simulation Monte Carlo (DSMC) technique, 

(b) . The Direct Solution of the Boltzmann Equation. 

The DSMC is a straight forward and intuitive method to acquire the flowfield in 
large Knudsen number (Bird, 1990). The DSMC method is capable of simulating the 
low-density flow phenomena involving nonequilibrium of the translational and internal 

modes. Calculations can be made throughout the rarefied regimes and can overlap with 

* 

the continuum calculations, the degree of overlap being dependent on the magnitude of 
the computating resources. It is very practical to adopt DSMC technique to attack the 
vacuum plume problems. 

Direct solutions of Boltzmann equation have been a difficult subject for the aerody- 
namic problems (Cercignani, 1988). Based on literature survey, the only two-dimensional 
flow simulations based on direct solutions of Boltzmann equation was published by 
Tcheremissine (1989). However, direct solutions of Boltzmann equation in multidimen- 
sional flows require great computer resources as well as new numerical techniques. We 
plan to solve the nonlinear Boltzmann equation in the high Knudsen number backflow 
region of the plume using ALE technique. The appropriate upstream boundary condi- 
tions (close to the nozzle exit with relatively larger Knudsen number) will be provided 
by the solution of Burnett equations obtained by DLUF method. The solution of Boltz- 
mann equation will provide indepth knowledge in classical sense of Hydrodynamics. 
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Derivation of the Burnett Equations in Cartesian Coordinates 


The Burnett equations in primitive variables are cast in the Cartesian tensor forms: 


+ P,iV i + — 0 

dVi 

dh dp 

P “b pVih } i P,i^i P Vi,i b - -f- 5i,i = 0 


(1) 

(2) 
(3) 


The constitutive relations for the viscous stress a t j, and heat flux q i%t , are 


_ _ ^,.(0) (1) (2) (*) 
a ij - °ij + a ij + <*ij +••• + %• 

(o) , (1) , (2) , , (k) 

Qi = q) + q) + q) + - • ■ + q) 


(4) 

(5) 


where i = 1,2,3 and j — 1,2,3, and the superscript number k = 0, l,2,...,n represent 
the Hh-order of accuracy. 

k=0, the Euler equations, 

k=l, the Navier-Stokes equations, 

k=2, the Burnett equations, 

k=3, the Super Burnett equations. 

The viscous stress and heat flux, up to second-order terms, are described as follows: 


dj 0) =P 8 ij (6) 

=-2p(D v -±D oa 6 l} ) (7) 

2 r\ -j -i 
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( 8 ) 


where 


-j~ OJ^R(A l j g ^-otoc &ij ) H“ ^4 ( B x j ^B aa 8ij) 

R 1 1 

^ ~C aa^ij) H" OJ^Eij ^ -^aor^iji ) } 


<T =o 

q ^ = -x™ 

q ' dxi 

( 2 ) _^ 2 so 1 dT du a 1 2 d du 


qr =-{0 


p XUl Tdx t dx a +02 T^3dx^ T dxj + 2 dx a dxj 

+ 6 *pdlt( Dai ~ 3 Dxx6at>> + 0A d^^ Dat ~ z Dxx8 ™) 

1 dT 1 

+ 0 5 -—( D at -- D xx 6 ai )}, 


_ 1 ( du x duj x 
ij ~ 2 { dx ] + dxj 

L . = i JL ( Aj?L) + 

,J 2 3xi pdxj 2 dxj pdxi 

M _ 1 <9« Q £?U a 

,} 2 dx a dxj dx Q dxi 

N'j = \[ ^f( D °i - - J^A.)] 


A„ = 


djT 

dx,dxj 


tJ 2 5xj 5xj dx x 
_ dT dT 
CiJ ~~d^ x Yx ~ 

T X j (^D xq , qT) wSiot^iD j a -D^^Sja). 


(9) 
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( 11 ) 

( 12 ) 

(13) 

(14) 

(15) 

(16) 

(17) 

(18) 
(19) 


where uq and 6 X can be determined by Chapman-Enskog method, depending on the gas 
molecular repulsive force models used. Expand the tensor terms from above equations 
in terms of the primitive variables and its derivatives, an explicit expression for the 
stress tensor components and heat flux vector components can be obtained. 
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